home *** CD-ROM | disk | FTP | other *** search
/ Celestin Apprentice 4 / Apprentice-Release4.iso / Languages / RLaB 1.18c / toolbox / bessely.r < prev    next >
Encoding:
Text File  |  1994-09-23  |  6.9 KB  |  249 lines  |  [TEXT/ttxt]
  1. ////-------------------------------------------------------------------//
  2.  
  3. //  Syntax:  bessely ( ALPHA, X )
  4.  
  5. //  Description:
  6.  
  7. //  Bessely is the bessel function of the second kind of order ALPHA.
  8. //  Either ALPHA or X may a vector in which case a vector result of
  9. //  the same dimension is returned.  If both ALPHA and X are vectors
  10. //  a matrix result of dimension length(ALPHA) x length(X) is returned.
  11. //  If either ALPHA or X is a scalar, the other argument may be a matrix
  12. //  and a matrix result of the same dimension is returned.
  13.  
  14. //-------------------------------------------------------------------//
  15.  
  16. // See Numerical Recipes in C (second edition)
  17. // Currently only integer order supported.
  18.  
  19. static (besselyn, bessely0, bessely1);
  20.  
  21. bessely = function(alph,x) 
  22. {
  23.   global (pi)
  24.  
  25.   // Checking the class will automatically result in an error if
  26.   // the argument doesn't exist.
  27.   if(class(alph) != "num" || class(x) != "num") {
  28.     error("bessely: argument is non-numeric.");
  29.   }
  30.  
  31.   if(min(size(alph)) > 1 && min(size(x)) > 1) {
  32.     // Both arguments are matrices.
  33.     error("bessely: only one argument may be a matrix.");
  34.   }
  35.  
  36.   if(length(alph) == 1) {
  37.     // x is a matrix (or scalar), alph a scalar
  38.     if(alph == 0 || mod(alph,int(alph)) == 0) {
  39.       return besselyn(alph,x);
  40.     else
  41.       error("bessely: fractional orders not yet supported.");
  42.       // return besselyr(alph,x);
  43.     }
  44.   else if(length(x) == 1) {
  45.     // alph is a matrix (or vector), x a scalar
  46.     y = zeros(size(alph));
  47.  
  48.     // Use integer order routines for integer alph
  49.     inte_ind = find(alph == 0 || mod(alph,int(alph)) == 0);
  50.     frac_ind = complement(inte_ind,1:alph.n);
  51.  
  52.     if(inte_ind.n) {
  53.       for(index in inte_ind) {
  54.         y[index] = besselyn(alph[index],x);
  55.       }
  56.     }
  57.     if(frac_ind.n) {
  58.       error("bessely: fractional orders not yet supported.");
  59.       // for(index in frac_ind) {
  60.         // y[index] = besselyr(alph[index],x);
  61.       // }
  62.     }
  63.     return y;
  64.   else
  65.     // alph and x are both vectors
  66.     y = zeros(length(alph),length(x));
  67.  
  68.     // Use integer order routines for integer alph
  69.     inte_ind = find(alph == 0 || mod(alph,int(alph)) == 0);
  70.     frac_ind = complement(inte_ind,1:alph.n);
  71.  
  72.     if(inte_ind.n) {
  73.       for(index in inte_ind) {
  74.         y[index;] = besselyn(alph[index],x);
  75.       }
  76.     }
  77.     if(frac_ind.n) {
  78.       error("bessely: fractional orders not yet supported.");
  79.       // for(index in frac_ind) {
  80.         // y[index;] = besselyr(alph[index],x);
  81.       // }
  82.     }
  83.     return y;
  84.   }}
  85. };
  86.  
  87. // In these functions x can be a vector (or matrix), but n
  88. // must be a scalar.
  89.  
  90. // No argument checking is performed on static functions.
  91.  
  92. besselyn = function ( n, x ) {
  93.  
  94.   local(abs_x,y,index,index_1,index_2,index_3,arg,p_1,p_2,p_3, ...
  95.         two_over_x,m,y_2,Norm,even,limit,lim_ind);
  96.  
  97.   if(n == 0) {
  98.     return bessely0(x);
  99.   else if (n == 1) {
  100.     return bessely1(x);
  101.   else
  102.     // integer order greater than two
  103.  
  104.     y = zeros(size(x));
  105.  
  106.     index_1 = 1:x.n;
  107.     index_2 = find(x == 0);
  108.  
  109.     index_1 = complement(intersection(index_1,index_2),index_1);
  110.  
  111.     if(index_2.n) {
  112.       y[index_2] = -1./zeros(size(index_2));
  113.     }
  114.  
  115.     if(index_1.n) {
  116.       arg = x[index_1];
  117.  
  118.       p_1 = bessely0(arg);
  119.       p_2 = bessely1(arg);
  120.       two_over_x = 2./arg;
  121.       for(index in 1:(n-1)) {
  122.         p_3 = index * (two_over_x .* p_2) - p_1;
  123.         p_1 = p_2;
  124.         p_2 = p_3;
  125.       }
  126.       y[index_1] = p_3;
  127.     }
  128.  
  129.     return y;
  130.   }}
  131. };
  132.  
  133. bessely0 = function(x) {
  134.   local(index_1,index_2,index_3,y,c_1,c_2,arg_1,arg_2,p_1,p_2);
  135.  
  136.   index_1 = find(x < 8);
  137.   index_2 = complement(index_1,1:x.n);
  138.   index_3 = find(x == 0.);
  139.  
  140.   index_1 = complement(intersection(index_1,index_3),index_1);
  141.  
  142.   // First evaluate using rational approx. for small arguments.
  143.  
  144.   y = zeros(size(x));
  145.  
  146.   if(index_3.n) {
  147.     y[index_3] = -1./zeros(size(index_3));
  148.   }
  149.  
  150.   if(index_1.n) {
  151.     c_1 = [ -2957821389, 7062834065, -512359803.6, ...
  152.             10879881.29, -86327.92757, 228.4622733 ];
  153.     c_2 = [ 40076544269, 745249964.8, 7189466.438, ...
  154.             47447.26470, 226.1030244 ];
  155.  
  156.     arg_1 = x[index_1];
  157.     arg_2 = arg_1 .* arg_1;
  158.  
  159.     // Rational approximation to the "regular part"
  160.     p_1 = c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + arg_2 .* ...
  161.           (c_1[4] + arg_2 .* (c_1[5] + arg_2 .* c_1[6]))));
  162.     p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* ...
  163.           (c_2[4] + arg_2 .* (c_2[5] + arg_2))));
  164.     y[index_1] = p_1./p_2 + (2/pi)*besselj(0,arg_1).*log(arg_1);
  165.   }
  166.  
  167.   // Now evaluate for large arguments using approximating form.
  168.  
  169.   if(index_2.n) {
  170.     c_1 = [ -0.1098628627e-2, 0.2734510407e-4, -0.2073370639e-5, ...
  171.             0.2093887211e-6 ];
  172.     c_2 = [ -0.1562499995e-1, 0.1430488765e-3, -0.6911147651e-5, ...
  173.             0.7621095161e-6, -0.934945152e-7 ];
  174.  
  175.     arg_1 = 8./x[index_2];
  176.     arg_2 = arg_1 .* arg_1;
  177.  
  178.     p_1 = 1 + arg_2 .* (c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + ...
  179.           arg_2 .* c_1[4])));
  180.     p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* (c_2[4] + ...
  181.           arg_2 .* c_2[5])));
  182.  
  183.     c_1 = x[index_2];
  184.     c_2 = c_1 - pi/4;
  185.     y[index_2] = sqrt(2./(pi*c_1)).*(p_1.*sin(c_2) + arg_1.*p_2.*cos(c_2));
  186.   }
  187.  
  188.   return y;
  189. };
  190.  
  191.  
  192. bessely1 = function(x) {
  193.   local(index_1,index_2,index_3,y,c_1,c_2,arg_1,arg_2,p_1,p_2);
  194.  
  195.   index_1 = find(x < 8);
  196.   index_2 = complement(index_1,1:x.n);
  197.   index_3 = find(x == 0.);
  198.  
  199.   index_1 = complement(intersection(index_1,index_3),index_1);
  200.  
  201.   // First evaluate using rational approx. for small arguments.
  202.  
  203.   y = zeros(size(x));
  204.  
  205.   if(index_3.n) {
  206.     y[index_3] = -1./zeros(size(index_3));
  207.   }
  208.  
  209.   if(index_1.n) {
  210.     c_1 = [ -0.4900604943e13, 0.1275274390e13, -0.5153438139e11, ...
  211.             0.7349264551e9, -0.4237922726e7, 0.8511937935e4 ];           
  212.     c_2 = [ 0.2499580570e14, 0.4244419664e12, 0.3733650367e10, ...
  213.             0.2245904002e8, 0.1020426050e6, 0.3549632885e3 ];
  214.     
  215.     arg_1 = x[index_1];
  216.     arg_2 = arg_1 .* arg_1;
  217.  
  218.     p_1 = arg_1.*(c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + arg_2 .* ...
  219.           (c_1[4] + arg_2 .* (c_1[5] + arg_2 .* c_1[6])))));
  220.     p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* ...
  221.           (c_2[4] + arg_2 .* (c_2[5] + arg_2 .* (c_2[6] +arg_2)))));
  222.     y[index_1] = p_1./p_2 + (2/pi)*(besselj(1,arg_1).*log(arg_1) - 1./arg_1);
  223.   }
  224.  
  225.   // Now evaluate for large arguments using approximating form.
  226.  
  227.   if(index_2.n) {
  228.     c_1 = [ 0.183105e-2, -0.3516396496e-4, 0.2457520174e-5, ...
  229.             -0.240337019e-6 ];
  230.     c_2 = [ 0.04687499995, -0.2002690873e-3, 0.8449199096e-5, ...
  231.             -0.88228987e-6, 0.105787412e-6 ];
  232.  
  233.     arg_1 = 8./x[index_2];
  234.     arg_2 = arg_1 .* arg_1;
  235.  
  236.     p_1 = 1 + arg_2 .* (c_1[1] + arg_2 .* (c_1[2] + arg_2 .* (c_1[3] + ...
  237.           arg_2 .* c_1[4])));
  238.     p_2 = c_2[1] + arg_2 .* (c_2[2] + arg_2 .* (c_2[3] + arg_2 .* (c_2[4] + ...
  239.           arg_2 .* c_2[5])));
  240.  
  241.     c_1 = x[index_2];
  242.     c_2 = c_1 - 3*pi/4;
  243.     y[index_2] = sqrt(2./(pi*c_1)).*(p_1.*sin(c_2) + arg_1.*p_2.*cos(c_2));
  244.  
  245.   }
  246.  
  247.   return y;
  248. };
  249.